home *** CD-ROM | disk | FTP | other *** search
/ Trading on the Edge / Trading On The Edge - CD-ROM Toolkit (Wayzata Technology)(2031)(1994).bin / pc / mac_file / software / nn_prepr / nlds.wgz < prev   
Text File  |  1992-09-21  |  15KB  |  488 lines
  1. Scripts for  Nonlinear Analysis of  Time-series
  2.  
  3. FUNCTION Hurst(closecolumn,rowstart,rowend,differences,N,
  4.         start_N,delta_N,rscolumn,report,graph)
  5. { script to calculate the Hurst coefficient and}
  6. {where}
  7. {    N  = the number of windows (number of graph points on (R/S,N)}
  8. {    N_start = initial value of N or size of first window}
  9. {    delta_N = arithmetic increment in window size}
  10. {        or if delta_N=0 then geometric increment in window size}
  11. {    Eu(t)  = the influx for each t}
  12. {    M(N)   = the average influx for each N}
  13. {    X(t,N) = the difference between Eu(t) and M(N)}
  14. {    R(N)   = the difference between the Max(Xt,N) and Min(Xt,N)}
  15. {    R/S or H = the rescaled range or R(N) / stdev.}
  16. {    b0        = the intercept of log(R/S) vs log(N)}
  17. {    b1        = the slope of log(R/S) vs log(N)}
  18. {    f        = fractal dimension of the time-series or 1/H}
  19. {    differences  = flag for calculating Eu over series or first}
  20. {        difference of series}
  21. {    report = flag for recording details of R/S computations}
  22. {    graph  = flag for plotting a scatter diagram of log(r/s) vs log(N)}
  23. {Hurst coefficient}
  24.     {-> if H >0.5 persistent trends,}
  25.     {-> if H = 0.5 random}
  26.     {-> if H <0.5 anti-persistent trends}
  27. Column Numbers
  28. Manual Recalc
  29. define i,j,k,l,counter,no_of_records,step_size,no_of_windows, rowcounter
  30. define M,S,Eu,Xt,R,RS,H, lnRS,lnN, f,b0,b1, Max_of_Xt, Min_of_Xt, totalXt
  31. define change_Eu,windowsize,deltaN, end_row
  32. no_of_records = rowend-rowstart+1
  33. step_size=start_N
  34. {column labels}
  35.     if report = 1
  36.     put "Eu" into makecell(rscolumn,1)
  37.     put "M" into makecell(rscolumn+1,1)
  38.     put "S" into makecell(rscolumn+2,1)
  39.     put "Xt" into makecell(rscolumn+3,1)
  40.     put "sumX" into makecell(rscolumn+4,1)
  41.     put "R" into makecell(rscolumn+5,1)
  42.     put "R/S" into makecell(rscolumn+6,1)
  43.     end if
  44. { calculate the differences from t to t+1}
  45. for i = rowstart+1 to rowend
  46.     if differences = 1
  47.     Eu = ln(indirect(makecell(closecolumn,i))
  48.         /indirect(makecell(closecolumn,i-1)))
  49.     put Eu into makecell(rscolumn,i)
  50.     else
  51.     put indirect(makecell(closecolumn,i)) into makecell(rscolumn,i)
  52.     end if
  53. end for
  54. FOR counter = 1 to N
  55.     {calculate for each N:
  56.     the number of windows,the windowsize starting with start_N}
  57.     if delta_N = 0
  58.     {if delta_N is put to zero then value of N will double}
  59.     windowsize = start_N^counter
  60.     else
  61.     {if start_n = 2 and delta_N = 2 then when N = 2,
  62.         windowsize = 2+ 2*(2-1) =4}
  63.     windowsize = start_N + delta_N * (counter-1)
  64.     end if
  65.     no_of_windows = int((no_of_records-1)/ windowsize)
  66.     put counter into makecell(rscolumn-1,1)
  67.     {calculate the last low for reach windowsize}
  68.     end_row = rowstart+(no_of_windows*windowsize)
  69.     FOR i = rowstart+1 TO  end_row step windowsize
  70.     {calculate M and S for the current windowsize}
  71.     put i into makecell(rscolumn-1,2)
  72.     M = avg(range(makerange(rscolumn,i, rscolumn,i+windowsize-1)))
  73.     S = stdev(range(makerange(rscolumn,i, rscolumn,i+windowsize-1)))
  74.     if report = 1
  75.         put M into makecell(rscolumn+1,i)
  76.         put S into makecell(rscolumn+2,i)
  77.     end if
  78.     max_of_Xt = -10000000
  79.     min_of_Xt =  10000000
  80.     totalXt= 0
  81.     {calculate for each N AND for each window Xt, R and RS}
  82.     FOR j = i TO i+windowsize-1
  83.         Eu = indirect(makecell(rscolumn,j))
  84.         {Calculate X(t,N)  or e(u)-M for N windows}
  85.         Xt = Eu - M
  86.         put Xt into makecell(rscolumn+3,j)
  87.         put j into makecell(rscolumn-1,3)
  88.         totalXt=totalXt + Xt
  89.         put totalXt into makecell(rscolumn+4,j)
  90.     END FOR { end of computations for a single window}
  91.     {calculate R (N) for N windows}
  92.     max_of_Xt=max((range(
  93.         makerange(rscolumn+4,i, rscolumn+4,i+windowsize-1))))
  94.     min_of_Xt=min((range(
  95.         makerange(rscolumn+4,i, rscolumn+4,i+windowsize-1))))
  96.     put max_of_Xt into makecell(rscolumn+7,i)
  97.     put min_of_Xt into makecell(rscolumn+8,i)
  98.     {calculate R, RS}
  99.     R = (max_of_Xt) - (min_of_Xt)
  100.     if S = 0
  101.         RS=R/0.000000000000001 {avoids division by 0}
  102.     else
  103.         RS = R/S
  104.     end if
  105.     put R into makecell(rscolumn+5,i)
  106.     put RS into makecell(rscolumn+6,i)
  107.     END FOR {end of computation for a single N}
  108.     {report summary results N, ln N, R/S, ln R/S}
  109.     put "data point" into makecell(rscolumn+10,1)
  110.     put " N " into makecell(rscolumn+11,1)
  111.     put " R/S " into makecell (rscolumn+12,1)
  112.     put " log R/S " into makecell (rscolumn+13,1)
  113.     put "log N" into makecell(rscolumn+14,1)
  114.     put "Column = " into makecell(rscolumn+15,1)
  115.     put closecolumn into makecell(rscolumn+16,1)
  116.     select range(makerange(rscolumn+10,1,rscolumn+16,1))
  117.     align center
  118.     text style "B"
  119.     unselect
  120.     put counter into makecell(rscolumn+10,counter+1)
  121.     { record N}
  122.     put windowsize into makecell(rscolumn+11,counter+1)
  123.     { record average or R/S values}
  124.     RS = avg(range(makerange(rscolumn+6,rowstart, rscolumn+6,rowend)))
  125.     put RS into makecell(rscolumn+12,counter+1)
  126.     put log(RS) into makecell(rscolumn+13,counter+1)
  127.     put log(windowsize) into makecell(rscolumn+14,counter+1)
  128.     { compute the H coefficient for each N}
  129.     {clean up the worksheet :erase sumX & RS values}
  130.     select range(makerange(rscolumn+1,rowstart, rscolumn+8,rowend))
  131.     clear
  132.     unselect
  133.     select range(makecell(rscolumn-1,1))
  134.     unselect
  135. END FOR {end of computations for all N's}
  136. { run regression of log(R/S) vs log (N)}
  137. counter = 0
  138. for i = 2 to N+1
  139.     if indirect(makecell(rscolumn+13,i)) <= 0.5
  140.     counter = counter +1
  141.     end if
  142. end for
  143. select range(makerange(rscolumn+13,2, rscolumn+14,counter))
  144. select more range(makerange(rscolumn+13,N+3, rscolumn+14,N+3))
  145. select more range(makerange(rscolumn+13,N+5, rscolumn+20,N+10))
  146. regress
  147. put "log(a)" into makecell(rscolumn+13,N+2)
  148. put "Hurst" into makecell(rscolumn+14,N+2)
  149. put "Correlation" into makecell(rscolumn+15,N+2)
  150. H = indirect(makecell(rscolumn+14,N+3))
  151. put 2^(2*H-1)-1 into makecell(rscolumn+15,N+3)
  152. select range(makerange(rscolumn+13,N+3, rscolumn+15,N+3))
  153. precision 4
  154. unselect
  155. {plot a scatter diagram of log(R/S) vs log(N)}
  156. if graph = 1
  157.     {copy the log R/S and log (N)}
  158.     select range(makerange(rscolumn+13,2,rscolumn+13,N+2))
  159.     copy
  160.     go to cell makecell(rscolumn+20,2)
  161.     paste
  162.     select range(makerange(rscolumn+14,2,rscolumn+14,N+2))
  163.     copy
  164.     go to cell makecell(rscolumn+19,2)
  165.     paste
  166.     {plot a scatter; adjust the scales on the axis}
  167.     Unselect
  168.     Add Chart Range frac(R1C25..R28C30,43,1,64,252)
  169.     Using range(makerange(rscolumn+19,2,rscolumn+20,N+2))
  170.     object name "R/Splot"
  171.     select object "R/Splot"
  172.     Scatter
  173.     select chart "R/Splot" axis 3
  174.     manual scaling from 0 to 1 with 10 major and 5 minor divisions
  175.     axis log base 10
  176.     select chart "R/Splot" axis 1
  177.     manual scaling from 0 to 5 with 5 major and 5 minor divisions
  178.     axis log base 10
  179.     Hide Legend
  180. end if
  181. { save file}
  182. save
  183. END FUNCTION
  184.  
  185.  
  186. FUNCTION CorrDim(col,rowstart,rowend,start_dimem,end_dimem,
  187.         tau,dt,R,rscolumn)
  188. { script to calculate the Correlation Dimension as}
  189. {     defined by Grassberger & Procaccia}
  190. { where}
  191. { npt        =    number of observations}
  192. { dimem        =   embedding dimensions}
  193. { tau        =    lag time for phase space}
  194. { dt        =    increase in each measurement}
  195. { R        =    the beginning distance}
  196. { CR        =    correlation integral}
  197. { lag        =    minimum time between pairs}
  198. { col        =    input data column}
  199. Column Numbers
  200. Manual Recalc
  201. Define i,j,k, npt, ind, total,its, z,d
  202. Define theta,theta2, CR, lag,rownumber, colnumber, L
  203. Define stop1,dimem
  204. Define loop
  205. {this loop provides for multiple estimates of corr dim}
  206. for loop = start_dimem to end_dimem
  207.     dimem = loop
  208.     if loop > start_dimem
  209.     rscolumn = rscolumn + 10
  210.     end if
  211.     theta = 0
  212.     theta2 = 0
  213.     CR = 0
  214.     ind = 1
  215.     k = 1
  216.     lag = 0
  217.     total = 0
  218.     its = 0
  219.     npt = rowend - rowstart + 1
  220.     rownumber =1
  221.     colnumber = 10
  222.     {read the input column and reconstruct the phase space}
  223.     for i = rowstart to rowstart+npt
  224.     for j = 1 to dimem
  225.         {read  x(i+(j-1) * tau and put it into a}
  226.         {  working area z(i,j) on the sheet}
  227.         z = indirect(makecell(col,i +(j-1)*tau))
  228.         put z into makecell(rscolumn+colnumber+j,i)
  229.     end for
  230.     end for
  231.     { calculate the maximum length of the phase space}
  232.     npt = npt -dimem*tau
  233.     stop1 = 0
  234.     while (stop1 = 0)
  235.     for k = rowstart to rowend
  236.         for i = rowstart to rowend
  237.         d = 0
  238.         for j = 1 to dimem
  239.             {calculating the square of distance}
  240.             d = d + (indirect(makecell(rscolumn+colnumber+j,i-lag))
  241.             - indirect(makecell(rscolumn+ colnumber+j,i))) ^2
  242.         end for
  243.         d = sqrt(d)
  244.         if d > R
  245.             theta2 = 0
  246.         else
  247.             theta2 = 1
  248.         end if
  249.         theta = theta + theta2
  250.         end for
  251.         lag = lag +1
  252.     end for
  253.     {calc correlation integral}
  254.     CR = (1/npt^2)*theta
  255.     {produce result table}
  256.     rownumber = rownumber + 1
  257.     put "Counter" into makecell(rscolumn,1)
  258.     put "Correlation Integral" into makecell(rscolumn+1 ,1)
  259.     put " Distance R" into makecell(rscolumn+2,1)
  260.     put "log(CR)" into makecell(rscolumn+3, 1)
  261.     put "log(R)" into makecell(rscolumn+4,1)
  262.     put rownumber-1 into makecell(rscolumn, rownumber)
  263.     put CR into makecell (rscolumn+1,rownumber)
  264.     put R into makecell (rscolumn+2, rownumber)
  265.     put log(cr) into makecell(rscolumn+3, rownumber)
  266.     put log(R) into makecell(rscolumn+4,rownumber)
  267.     L = L + 1
  268.     if L > 20
  269.         stop1 =1
  270.     else
  271.         R = R + dt
  272.         CR = 0
  273.         theta = 0
  274.         theta2 = 0
  275.         lag = 0
  276.     end if
  277.     end while
  278.     {format results}
  279.     select range(makerange(rscolumn, 1, rscolumn,rownumber))
  280.     precision 0
  281.     select range(makerange(rscolumn+1, 1, rscolumn+4,rownumber))
  282.     precision 6
  283.     unselect
  284.     { clean up spreadsheet}
  285.     select range(makerange(
  286.     rscolumn+colnumber, rowstart, rscolumn+colnumber+dimem,rowend))
  287.     clear
  288.     unselect
  289.     { calc the regression}
  290.     select range(makerange(rscolumn+3, 2, rscolumn+4,rownumber))
  291.     select more range(makerange(
  292.     rscolumn+3, rownumber+2, rscolumn+4,rownumber+2))
  293.     select more range(makerange(
  294.     rscolumn+3, rownumber+4, rscolumn+10,rownumber+10))
  295.     regress
  296.     select range(makerange(
  297.     rscolumn+3, rownumber+2, rscolumn+5,rownumber+12))
  298.     precision 6
  299.     unselect
  300. end for
  301. END FUNCTION
  302.  
  303.  
  304. FUNCTION Lyapunov(col,rowstart,rowend,start_dimem,end_dimem,tau,dt,
  305.             scalmax,scalmin,evolve,lag,rscolumn)
  306. { script to calculate the Lynapunov exponent    }
  307. { where}
  308. { npt        =    number of observations}
  309. { dimem        =   embedding dimensions}
  310. { tau        =    lag time for phase space}
  311. { dt        =    increase in each measurement}
  312. { scalmax   =    maximum divergence}
  313. { scalmin   =    minimum divergence}
  314. { evolve    =    evolution time}
  315. *{ lag        =    minimum time between pairs}
  316. { col        =    input data column}
  317. { wrkcol    =    first working column}
  318. { wrkrow    =    first working row cell}
  319. Column Numbers
  320. Manual Recalc
  321. Define i,j, npt, ind,ind2,indold, total,its, pt1,pt2, z
  322. Define zlyap, zmult, anglmx,thmin,iii,dnew,dot,cth,d,di,dii,df,th
  323. Define dummy,rownumber, dimem
  324. Define stop1,stop2,loop
  325. for loop = start_dimem to end_dimem
  326.     dimem = loop
  327.     if loop > start_dimem
  328.     rscolumn = rscolumn +1
  329.     end if
  330.     npt = rowend - rowstart +1
  331.     ind = 1
  332.     total = 0
  333.     its = 0
  334.     rownumber = 1
  335.     {read the input column and reconstruct the phase space}
  336.     for i = rowstart to rowstart + npt-(dimem-1)*tau
  337.     for j = 1 to dimem
  338.         {read  x(i+(j-1) * tau and put it into a
  339.         working area on the sheet}
  340.         z = indirect(makecell(col,i +(j-1)*tau))
  341.         put z into makecell(rscolumn+10+j,i)
  342.     end for
  343.     end for
  344.     {max length of phase space}
  345.     npt = npt -dimem*tau-evolve
  346.     di = 1000000000
  347.     {find initial pair}
  348.     for i = rowstart+lag +1 to rowstart+npt
  349.     d = 0
  350.     for j = 1 to dimem
  351.         {calc distance}
  352.         d = d + (indirect(makecell(rscolumn+10+ind,i))
  353.           - indirect(makecell(rscolumn+10+j,i)) )^2
  354.     end for
  355.     d =sqrt(D)
  356.     if d > di  or d < scalmin
  357.         continue for { store the best point}
  358.     end if
  359.     di = d
  360.     ind2 = i
  361.     end for
  362.     {coordinates of evolve}
  363.     stop1 = 0
  364.     stop2 = 0
  365.     while (stop1 = 0)
  366.     for j = 1 to dimem
  367.         pt1 = indirect(makecell(j,ind+evolve))
  368.         put pt1 into makecell(rscolumn+20,j)
  369.         pt2 = indirect(makecell(j, ind2+evolve))
  370.         put pt2 into makecell(rscolumn+21,j)
  371.     end for
  372.     df = 0
  373.     {compute final divergence}
  374.     for j = 1 to dimem
  375.         df = df+ (indirect(makecell(rscolumn+21,j))
  376.         - indirect(makecell(rscolumn+20,j)) )^2
  377.     end for
  378.     df = sqrt(df)
  379.     its = its +1
  380.     total = total + (log(df/di)/(evolve*dt*log(2)))
  381.     zlyap =total/its
  382.     {produce output table}
  383.     rownumber = rownumber +1
  384.     put "zlyap" into makecell(rscolumn,1 )
  385.     put "evolve*its" into  makecell(rscolumn+1,1 )
  386.     put "di" into  makecell(rscolumn+2,1 )
  387.     put "df" into  makecell(rscolumn+3,1 )
  388.     put zlyap into makecell(rscolumn,rownumber)
  389.     put evolve*its into  makecell(rscolumn+1,rownumber )
  390.     put di into  makecell(rscolumn+2,rownumber )
  391.     put df into  makecell(rscolumn+3,rownumber )
  392.     indold = ind2
  393.     zmult = 1
  394.     anglmx =.3
  395.     stop2=0
  396.     while (stop1 = 0) and (stop2 = 0)
  397.         thmin = 3.14
  398.         {look for replacement points}
  399.         for i =rowstart to rowstart+npt
  400.         iii = abs(int(i-(ind+evolve)))
  401.         {reject if replacement point is too close to original}
  402.         if iii < lag
  403.             continue for
  404.         end if
  405.         dnew = 0
  406.         for j = 1 to dimem
  407.             dnew = dnew+ (indirect(makecell(rscolumn+20,j))
  408.             - indirect(makecell(rscolumn+10+j,i)) ) ^2
  409.         end for
  410.         dnew = sqrt(dnew)
  411.         if (dnew > zmult * scalmax) or ( dnew < scalmin)
  412.             continue for
  413.         end if
  414.         dot = 0
  415.         for j = 1 to dimem
  416.             {calc pt1 - pt2}
  417.             dummy = indirect(makecell(rscolumn+20,j))
  418.               -indirect(makecell(rscolumn+21,j))
  419.             dot = dot + (indirect(makecell(rscolumn+20,j))
  420.               -indirect(makecell(rscolumn+10+j,i)) * dummy)
  421.         end for
  422.         cth = abs(dot/dnew*df)
  423.         if cth > 1
  424.             cth =1
  425.         end if
  426.         th = cos(cth)
  427.         if th > thmin
  428.             continue for
  429.         end if
  430.         thmin = th
  431.         dii = dnew
  432.         ind2 = i
  433.         end for
  434.         if (thmin < anglmx)
  435.         ind =ind + evolve
  436.         if (ind >= npt)
  437.             stop1 = 1
  438.             stop2 = 1
  439.         else
  440.             di=dii
  441.             stop2 = 1
  442.             continue while
  443.         end if
  444.         else        {case when thmin => anglmx}
  445.         zmult = zmult + 1
  446.         if (zmult < 5)
  447.             stop2 = 0
  448.             continue while
  449.         else  { case when zmult => 5}
  450.             zmult = 1
  451.             anglmx = 2*anglmx
  452.             if (anglmx < 3.14)
  453.             stop2 = 0
  454.             continue while
  455.             else     {case when anglmx => 3.14}
  456.             ind2=indold+evolve
  457.             dii=df
  458.             ind=ind+evolve
  459.             if ( ind >= npt)
  460.                 stop1 = 1
  461.                 stop2 = 1
  462.             else
  463.                 di = dii
  464.                 stop2 = 1
  465.                 continue while
  466.             end if
  467.             end if
  468.         end if
  469.         end if
  470.     end while
  471.     end while
  472.     {clean up the worksheet}
  473.     select range(makerange(rscolumn+10,rowstart,rscolumn+21,rowend))
  474.     clear
  475.     unselect
  476.     {run regression on lynapunov exponent}
  477.     select range(makerange(
  478.         rscolumn,rownumber-20,rscolumn+1,rownumber))
  479.     select more range(makerange(
  480.         rscolumn,rownumber+2,rscolumn+1,rownumber+2))
  481.     select more range(makerange(
  482.         rscolumn,rownumber+4,rscolumn+5,rownumber+5))
  483.     regress
  484. end for {for loop}
  485. Return("done")
  486. END FUNCTION
  487.  
  488.